.. _tutorial_6:
TUTORIAL 6 : Free Energy (FE) and replica-exchange methods
==========================================================
:Author: Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)
General concept
---------------
Most often it is impossible to calculate the absolute free energy value (we will see why later on).
However, *relative* free energy can always be expressed via the probability of the associated state:
.. math::
{\rm Free\ energy\ of\ state}\ = \ -kT \ln{\rm (probability\ of\ state)} \ \ \{+\ unknown\ const\}
then for the **FE difference** we get
.. math::
{\rm FED}_{12} = -kT \ln\frac{\rm probability\ of\ state\ 2}{\rm probability\ of\ state\ 1}
This is the most useful statistical-mechanical approach to FE evaluation, and all *advanced* FE methods are based on the above formula (althogh there exist also *thermodynamic* approaches, usually more demanding and time-consuming).
Example: Potential of mean force as FED
----------------------------------------
Potential of mean force (PMF or POMF) is defined as the work done upon reversibly bringing two particles from the state of zero-interaction into the state of full-interaction at a given distance, *W(R)*.
PMF can be directly obtained from RDF:
.. math::
\beta W(R) = -\ln g(R) + \ln g(\infty)
where the zero-interaction state corresponds to infinite separation :math:`(R\to\infty)`.
Alternatively, one can calculte
.. math::
\beta \Delta W(R_1 \to R_2) = -\ln g(R_2) + \ln g(R_1) = -\ln \frac{g(R_2)}{g(R_1)}
- the work to bring two particles :math:`(i,j)` from :math:`|\Delta\mathbf r_{ij}|=R_1` to :math:`R_2`.
Clearly, if :math:`g(R_1)=1`, we get back the total work :math:`W(R=\infty\to R_2)`.
.. figure:: Images-RDF/Ex2-gnuplot-PMF-prod1.png
:width: 640px
Why PMF is a free energy?
^^^^^^^^^^^^^^^^^^^^^^^^^
First, note that
.. math::
\frac{g(R_2)}{g(R_1)} =
\frac{ \langle \delta(|\Delta\mathbf r_{ij}|-R_2) \rangle }
{ \langle \delta(|\Delta\mathbf r_{ij}|-R_1) \rangle } =
\frac{ \int \delta(|\Delta\mathbf r_{ij}|-R_2)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} }
{ \int \delta(|\Delta\mathbf r_{ij}|-R_1)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} }
We can rewrite the canonical configuration integral (configurational part of the partition function)
.. math::
Q(N,V,T) = \int dR \int \delta(|\Delta\mathbf r|-R)\exp[-\beta U(\mathbf r^{N})] \ d{\mathbf r}^{N} =
\int Q_{NVT}(R) dR
.. Q(N,V,T) = \int \exp[-\beta U(\mathbf r^{N})] \ d\mathbf r^{N} =
as an integral (i.e. sum) over the substates corresponding to :math:`|\Delta\mathbf r_{ij}|=R`.
This implies that we can consider each substate as an individual thermodynamic state, whereas :math:`R`
serves as an additional **external constraint**, similar to *N,V,T* or *p*. Thereby,
.. math::
\beta \Delta W(R_1 \to R_2) = -\ln \frac{g(R_2)}{g(R_1)} =
\frac{ \langle \delta(|\Delta\mathbf r_{ij}|-R_2) \rangle }
{ \langle \delta(|\Delta\mathbf r_{ij}|-R_1) \rangle } = -\ln \frac{Q_{NVT}(R_2)}{Q_{NVT}(R_1)} = -\ln \frac{p(R_2)}{p(R_1)}
where
.. math::
p(R) = \frac{Q_{NVT}(R)}{Q(N,V,T)}
Therefore, we can evaluate the FED between any two substates (different particle separations) by performing a Metropolis
sampling along *R* and collecting a histogram of hits (or visits) on a discrete grid with sufficiently small bin size,
:math:`\Delta R`, which is of course equivalent to calculating :math:`g(R)`!
NOTE: One does not need to know the total partition function (nor the total configuration integral)
as it only serves as a normalization constant which cancels out!
- Question: what is the *"unknown constant"* in the first equation above?
- Similar formulae hold true generally, for other constrained substates (see the Appendix below).
FED methods in DL_MONTE
-----------------------
In general, Metropolis sampling is prone to getting trapped (stuck) in local FE minima and, hence,
suffering from undersampling (even avoiding!) regions around energetic or entropic barriers. This makes estimating FED
via probablity ratio troublesome...
The solution to this problem is to **bias the sampling** away from the minima and towards the barriers.
.. figure:: Images-FED/Slide2-FED1pic2.png
:width: 400px
.. figure:: Images-FED/Slide2-FED1.png
.. figure:: Images-FED/Slide2-FED2pic12a.png
:width: 800px
- **Ideally the bias potential is to perfectly compensate for the underlying free energy differences!**
- **This implies that the resulting probablity distribution in a perfectly biased simulation is uniform, i.e. flat!**
Parallel tempering with replica-exchange (RE)
---------------------------------------------
**Two approaches to parallelization in DL_MONTE**
.. figure:: Images-FED/Slide1-PAR-RE.png
In a **parallel multicanonical simualtion** each workgroup treats the system at a different tempreature
and generates it own MC trajectory. It is then possible to arrange periodic exchnage of configurations between
the workgroups, which is equivalent to concurrent jumps in *T* with the acceptance probability:
.. math::
acc[ \{ \mathbf{r}^N; T_1 \} \leftrightarrow \{ \tilde{\mathbf r}^N; T_2 \} ] =
\min(1, \exp \{- (\beta_1 U_{2} - \beta_2 U_2) - (\beta_2 U_1 - \beta_1 U_1) ] \} ) =
\min(1, \exp \{- \Delta \beta_{12} \Delta U_{12} \} )
If a FED evaluation is carried out in the same simulation, say FED(R), the acceptance rule includes an extra term:
.. math::
acc[ \{ \mathbf{r}^N; T_1 \} \leftrightarrow \{ \tilde{\mathbf r}^N; T_2 \} ] =
\min(1, \exp \{- \Delta \beta_{12} \Delta U_{12} - \Delta \beta_{12} \Delta U_{12}^{(bias)} \} )
Exercises
---------
Now try the exercises to learn how to perform FED calculations with DL_MONTE:
:ref:`tut6_ext1` - Reproducing interaction potentials *U(r)* with EE & WL/RE methods
:ref:`tut6_ext2` - Calculating PMF for two charged colloids in ionic cloud
:ref:`tut6_ext3` - Umbrella sampling in subranges (windows) with the use of WHAM method
Appendix: FED via a **Generalised (Expanded) Ensemble**
-------------------------------------------------------
The biasing is best illustrated with a *generalised, so-called Expanded Ensemble (EE) approach*.
The EE is defined by a generalised partition sum:
.. math::
Z_{EE} = \int w({\mathbf q}) d{\mathbf q} \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q})] d{\mathbf r}^N
or in a discretised form:
.. math::
Z_{EE} = \sum_m w({\mathbf q_m}) \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q_m})] d{\mathbf r}^N
where :math:`\{{\mathbf q}\}` stand for one or more **generalised coordinates (or constraints)**, which can also include normal thermodynamic variables (e.g. *N,V,T, p* etc), and :math:`w({\mathbf q})` is a **biasing weight** assigned to a particular
set of :math:`{\mathbf q}`.
Clearly, if Metropolis sampling is carried out in such a generalised (so-called expanded) ensemble, including MC moves
in :math:`{\mathbf q}`, the probablity of an EE substate
.. math::
p({\mathbf q}) = \frac{w({\mathbf q}) d{\mathbf q} \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q})] d{\mathbf r}^N}{Z_{EE}} =
\frac{w({\mathbf q}_m) \int_{ens} \exp[-\beta H({\mathbf r}^N; {\mathbf q}_m)] d{\mathbf r}^N}{Z_{EE}}
and the *effective* (biased) free energy difference
.. math::
\beta \Delta F_{nm}^{(EE)} = -\ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)} =
-\ln \frac{w({\mathbf q}_m) Q_{ens}({\mathbf q}_m)}{w({\mathbf q}_n) Q_{ens}({\mathbf q}_n)} =
-\ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} + \beta \Delta F_{nm}^{(true)}
Finally, if the weights :math:`\{w({\mathbf q}_m)\}` are chosen so that :math:`p({\mathbf q}_m)/p({\mathbf q}_n) = 1`,
we get:
.. math::
\beta\Delta F_{nm}^{(true)} - \ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} = 0
.. .. math::
.. \beta\Delta F_{nm}^{(EE)} = 0 \ \ \ {\rm and} \ \ \ \beta\Delta F_{nm}^{(true)} = \ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)}
- **Ideally the EE weights are to perfectly compensate for the unbiased probablity distribution!**
In practice, though, one has to arrange an iteration to self-consistently update the weights until
all the substates are populated reasonably well, after which the FED is obtained as
.. math::
\beta\Delta F_{nm}^{(true)} = \ln \frac{w({\mathbf q}_m)}{w({\mathbf q}_n)} - \ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)} =
\ln \frac{w({\mathbf q}_m) p({\mathbf q}_n)}{w({\mathbf q}_n)p({\mathbf q}_m)}
NOTE: It is convient (and conventional) to define the EE weights via a *biasing potential*:
.. math::
w({\mathbf q}) = \exp[-\beta U^{(bias)}({\mathbf q})]
so that
.. math::
\beta\Delta F_{nm} = -\Delta U_{nm}^{(bias)}({\mathbf q}) - \ln \frac{p({\mathbf q}_m)}{p({\mathbf q}_n)}
- **Ideally the bias potential is to perfectly compensate for the underlying free energy differences!**